function[L] = season_filter(fday,st_m,st_d,en_m,en_d),
%Filters fday (in mjd wrt 2000) for seasons
%st_m start month
%st_d = start day
%en_m = end month
%en_d = end day
%Returns Logical array L
%Limitations !!
% careful while giving Feb 28, Feb29 etc. No check is done for leap year
% Suggestion for feb last days is just use March 1st !!
%Also seasons that spans over the year boundaries (winter)
%Tested OK 05-Oct-2007


mjd_fday = datenum(2000,1,1); %find out the min and max year
c= datevec(min(fday)+mjd_fday);
minyear = c(1);

c= datevec(max(fday)+mjd_fday);
maxyear = c(1);

L = zeros(size(fday));

if st_m <= en_m,
    
    for i = minyear:maxyear,

    L = L| (fday >= datenum(i,st_m,st_d)-mjd_fday & fday <= datenum(i,en_m,en_d)-mjd_fday);

    end;
else, %if start month is less than end month (winter)
    
    for i = minyear:maxyear,

    L = L | (fday >= datenum(i,st_m,st_d)-mjd_fday & fday <= datenum(i,12,31)-mjd_fday);
    L = L | (fday >= datenum(i,1,1)-mjd_fday & fday <= datenum(i,en_m,en_d)-mjd_fday);

    end;
end;